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Abstract. In the near future the energy region above 
few hundreds of TeV may really be accessible for 
measurements of the atmospheric muon spectrum 
by the IceCube array. Therefore one expects that 
muon flux uncertainties above 50 TeV, related to a 
poor knowledge of charm production cross sections 
and insufficiently examined primary spectra and 
composition, will be diminished. We give predictions 
for the very high-energy muon spectrum at sea level, 
obtained with the three hadronic interaction models, 
taking into account also the muon contribution due 
to decays of the charmed hadrons. 
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I. Introduction 

The atmospheric neutrino flux at high energies is 
inevitably dominated by the prompt component due to 
decays of the charmed hadrons (D^, D°, D°, Df 
A+, . . .), hence the prompt neutrino flux becomes the 
major source of the background in searches for a diffuse 
astrophysical neutrino flux HI, Q, Q, (3- Insufficiently 
explored processes of the charm production give rise 
to most uncertainty in muon and neutrino fluxes. Ice- 
Cube, the first to begin operating as the Km3 neutrino 
telescope, has the real capability J5) to measure the 
atmospheric muon spectrum at energies up to 1 PeV and 
to shed light on the feasible range of the cross sections 
for the charmed particle production. 

In this work we try to extend to higher energies the 
muon flux calculations basing on the active hadronic in- 
teraction models with use of using the most reliable data 
of the primary cosmic ray measurements. We present 
results of the conventional muon flux calculations in the 
range 10-10 7 GeV using hadronic models QGSJET-II 
03 0, SIBYLL 2.1 Q, EPOS as well as Kimel 
and Mokhov (KM) [9], that were tested also in recent 
atmospheric muon flux calculations ifTOl . ifTTI . At the 
muon energies > 50 TeV we add the prompt muon 
flux originating from decays of the charmed hadrons, 
produced in collisions of cosmic rays with nuclei of air 
(for review see e.g. 02, 02, 03, 02). 



II. The method 

The high-energy muon fluxes are calculated using the 
approach |[T6ll to solve the atmospheric hadron cascade 
equations taking into account non-scaling behavior of 
inclusive particle production cross-sections, rise of to- 
tal inelastic hadron-nuclei cross-sections, and the non- 
power law primary spectrum (see also ifTTTn . 

To obtain the differential energy spectra of protons 
p(E, h) and neutrons n(E, h) at the depth h one needs 
to solve the set of equations: 

dN ± (E,h) N^E,^ 



dh 



Xn(E) J 



<P± N (E,x)N ± (E/x,h) ,, 



dx 
x 2 
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where N±(E, h) = p(E, h) ± n(E, h), 
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d<7pp(E , E) da pn (E , E) 



dE 
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dE 



\n{E) = 1/ [Noa l I^(E)j is the nucleon interaction 
length; x = E/Eq is the fraction of energy carried 
away by the secondary nucleon; da a b /dE - differential 
cross sections for inclusive reaction a + A — » b + X. 
The boundary conditions for Eq. ([TJ are N ± (E,0) — 
p {E)±n a (E). 

Consider that the solution of the system is 



iV ± (£', h) = 7V ± (£;,0)exp 



h(l-Z± N {E,h)) 



X N {E) 



(2) 



where Z^ N (E, h) are unknown functions. Substituting 
Eq. (f2]i into Eq. (Q]) we find the equation for these 
functions Z^ N (Z-factors): 



9{hZ% N ) 
dh 



^ N {E,x)^ N {E,x)x 



x exp [— tD NN [E, x, £)] dx, 
where rfc N (E,x) = x- 2 N ± {E/x,0)/N ± (E,0), 
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By integrating Eq. ^ we obtain the nonlinear integral 
equation 



dt 



dx$% N {E,x)r)% N {E,x)x 

xexp[-tD% N (E,x,t)], (4) 

which can be solved by iterations. The simple choice 
of zero-order approximation is Z^ffl (E, h) — 0, that is 

Dp$(E,x,h) = l/\ N (E/x)-l/\ N (E). For the n-th 
step we find 



1 

h 

x exp 



D%^(E,x,h) = 



h pi 

dt I dx$% N (E,x)i]± N (E,x)x 
Jo 

(5) 



(6) 

l-Zffi(E/x,h) l-Zp^(E,h) 
X N (E/x) X N (E) 

Nontrivial structure of the nucleon Z-factors (Fig. [U 
results from the non-power law behavior of the ATIC- 
2 primary spectrum (see later), non-scaling behavior of 
the particle production cross-sections and the energy 
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Fig. 1: The energy dependence of Z^ N {E, /i)-factors 
at the depth of 200 g cm~ 2 calculated for the ATIC-2 
primary spectrum. 



dependence of inelastic hadron-nucleus cross-sections. 
After that, using the obtained nucleon fluxes, we are 
able to calculate successively the meson and the lepton 
fluxes ifTTl . 

In our calculations we rely on recent data on the 
primary cosmic ray (PCR) spectra and composition 
obtained with Advanced Thin Ionization Calorimeter, 
balloon-borne experiment ATIC-2 IfTTl . In order to ex- 
tend the calculations to higher energies, up to 10 PeV, 
we use the data of the GAMMA experiment £18). 
The energy spectra and elemental composition, obtained 
in the GAMMA experiment, cover the 10 3 -10 5 TeV 
range and agree with the corresponding extrapolations 
of known balloon and satellite data at E > 10 3 TeV. 
Alternative primary spectra, used in the calculations for 
a region of very high energies, is the model by Zatsepin 
and Sokolskaya (ZS) [19|. The ZS proton spectrum at 
E > 10 6 GeV is compatible with KASCADE data J20) 
and the helium one is within the range of the KASCADE 
spectrum recontructed with use of QGSJET 01 and 
SIBYLL models. Besides we use also the Gaisser and 
Honda spectra (GH) [|2T| to compute the muon flux in 
more narrow energy range E^ < 10 5 GeV. 

III. High-energy muon spectra 

Apart from evident sources of AM, ir^ and K^a de- 
cays, we take into consideration three-particle semilep- 
tonic decays, Kr^, K® 3 . Also we take into account small 
fraction of the muon flux originated from decay chains 

K -> 7T -> H (K° s -> 7T+ + 7T - , -> 7T ± + 7T°). 



The prompt muon contribution due to decays of 
charmed hadrons at high energies is considered for 
the three charm production models, the recombination 
quark-parton model (RQPM) lfl2l . quark-gluon string 
model QGSM fl22), CTD, E), and the model by 
Pasquali, Reno and Sarcevic (PRS) [23] (see also 03)). 
Both RQPM and QGSM are the nonperturbative models, 
whereas PRS is based on the next-to-leading-order QCD 
calculations. 

The high energy spectra of conventional and prompt 
muons at ground level calculated for the vertical direc- 
tion are shown in Fig. |2] along with experimental data. 
The shaded areas here indicate the calculations of the 
conventional muon flux with KM model for the case of 
the ATIC-2 primary spectrum (light band on the left) and 
the GAMMA one (dark band on the right). The size of 
the bands correspond to statistical errors in the ATIC-2 
and GAMMA experiments. Solid, dashed, dash-dotted, 
and dotted lines indicate the calculations with use of 
ZS spectrum and set of hadronic models, KM, SIBYLL 
2.1, EPOS 1.61, and QGSJET-II. The experimental data 
comprise the measurements of L3+Cosmic 11241 . Cosmo- 
ALEPH ll25l as well as the data (converted to the 
surface) of deep underground experiments MSU (26), 
MACRO (27), LVD (28), Frejus (29), Baksan QUI . 
Artyomovsk ll3D . Notice that the calculation results do 
not fit well the Frejus and MSU data even though taking 
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Fig. 2: The high energy vertical muon spectra at ground level. The dashed-line curves and shaded areas present this 
work calculations with the ATIC-2 primary spectrum HT) < 10 TeV) and GAMMA one (E^ > 10 TeV). The 
solid curve mark the computation for the primary spectrum model by Zatsepin and Sokolskaya [19]. The numbers 
near thin lines indicate the prompt muon flux model: 1 - RQPM [12], 2 and 3 - the two calculation versions of 
PRS [23], 4 - QGSM fl2, lfl3l . 



into account the prompt muon component (thin lines in 
Fig. |2j. While, the LVD data are well described. 

The ZS model seems to be a reasonable bridge from 
TeV energy range to PeV one, providing a junction of 
the different energy ranges. Above 10 6 GeV the muon 
flux is apparently affected by the primary cosmic ray 
ambiguity in the vicinity of 'knee'. 

IV. Muon charge ratio 

The muon charge ratio depends on the proton to 
neutron ratio in primary cosmic rays as well as on the 
hadron production cross-sections. Thus, a comparison of 
the calculated ratio with experimental data in a 

wide energy range gives a possibility to study indirectly 
these features. At present time the muon charge ratio 
is measured with new facilities, that supply with high 
quality data from large number of muon events at high 
energies. 

In Fig. [3] we present our calculations of fi + / yT ratio 
along with the data of experiments [24], [25], l32l . 

E), El, El, ED, ED, IH, ED, Ho), ED. The 

calculations are made with the three hadronic interaction 
models and the two primary cosmic ray spectra. Solid 
line and dashed one mark the KM + GH computation 
for the zenith angles 0° and 90°, thin line indicates the 
KM + ZS result at 0°. Bold dotted line indicates the 



SIBYLL 2.1 + GH result, dotted shows the QGSJET- 
II + GH one, both for 0°. One may see that the 
KM and SIBYLL calculations reproduce data well with 
both versions of the primary spectra, GH and ZS. The 
calculated curves correspond approximately to the value 
1.3, that is in agreement with the recent data of BESS- 
TeV, CosmoALEPH, and L3+Cosmic experiments up to 
1 TeV. For higher energies the kaon source of muons 
becomes more intensive leading to a maximum of the 
muon ratio, ~ 1.4, at energy close to 10 TeV. This 
value agrees with the recent results of the MINOS far 
detector 11391 . The calculations also are in agreement 
with the spectrograph MUTRON data (35) at 9 = 89° 
including the point above 10 TeV. 

The QGSJET-II model (lower line) shows visible 
deviation from others: the predicted ratio is 

close to ~ 1.2, that might be explained by the higher 
extent of the proton and neutron flux equalization in the 
atmosphere due to reactions pA — > nX in the model 
(see Fig. [TJ. 

V. Conclusions 

One may expect that the primary cosmic ray uncer- 
tainty beyond the 'knee' would be negligible, if real 
prompt muons at energies around ~ 10 6 GeV (Fig. |2]i 
would appear at the flux level that not to be much lower 
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Fig. 3: Muon charge ratio at ground level computed for the three hadronic interaction models and the two primary 
cosmic ray spectra. Solid line marks the KM + GH result for 9 = 0°, dashed line shows the same at 90°. Thin 
line: the KM + ZS at 0°, bold-dotted: the SIBYLL 2.1 + GH, dotted (the lower): the QGSJET-II + GH at 0°. 



than the QGSM prediction. 

It seems reasonable to consider that the prompt muon 
flux higher than it is predicted by RQPM, is ex- 
cluded [3|. One may expect that more strong restriction 
of the prompt muon flux range will be extracted from 
the experiment in the near future. 
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